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Abstract 



This paper studies a two-dimensional chemotactic model for two species in 
which one of them produces a chemo-repellent for the other. It is shown asymp- 
totically and numerically how the chemical inhibits the invasion of a moving 
front for the second species and how stable steady states, which depend on the 
chemical concentration, can be reached. The results qualitatively explain exper- 
imental observations by Swain and Ray (Microbiol. Res. 164(2), 2009), where 
colonies of bacteria produce metabolite agents which prevent the invasion of 
fungi. 

Keywords: Chemotaxis; biocontrol; nonlinear dynamics; metastable patterns; 
front propagation. 

1 Introduction 

The movement of biological cells or microscopic organisms in response to chemical 
gradients is known as chemotaxis [2j [31] . It pertains to many biological phenomena, 
such as the motility of bacteria toward certain chemicals [39] , the response of fungal 
zoospores in the presence of metabolites [22 j, or the movement of endothelial cells 
toward angiogenic factors released by tumors [12 , among others. Since the semi- 
nal work by Keller and Segel [24[ [25 , mathematical modeling through systems of 
partial differential equations has established itself as an efficient tool to describe the 
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macroscopic self-organization patterns of cells occurring in chemotaxis (see [19] for a 
review) . 

This paper studies a chemotactic mathematical model describing the dynamics of 
two species (or colonies) of cells in which one of them releases a chemo-repellent for the 
other. Its design and study were motivated by the experimental results of Swain and 
Ray [37], who reported the beneficial activities (known as biocontrol) of the bacteria 
Bacillus subtilis in the presence of several harmful microflora found in cowdung, such 
as the phytopathogenic fungus Fusarium oxysporum. In one of the experiments, 
a uniform concentration of the F. oxysporum in a Petri dish was inoculated with 
the bacteria in two symmetrically located points along the diameter of the dish. 
The authors observed the inhibition of the in vitro growth of the fungus and the 
emergence of isolated patterns or strains, free of F. oxysporum concentration, which 
were localized near the places where the B. subtilis was applied. These patterns had 
the form of two circular fronts (see, for example, Figure 2 in [37]). The bacteria 
suppressed the growth of the fungus only after 48 hrs. of incubation and the isolation 
occurred in a stable manner for more than six days (144 hrs.); in other words, the 
fronts got stabilized and persisted after a long time, suggesting the emergence of 
stable transient states. 

It is a well-known fact [28] [32] that the B. subtilis produces antifungal metabolites 
(such as mycosubtilin) endowed with biocontrol properties against fungi like the F. 
oxysporum (see, for example, Leclere et al. [27]). The conditions under which the 
bacteria and the fungus isolate each other are, however, far from being understood. 
In order to describe the dynamics of the triggering of fungal suppression by the pres- 
ence of the metabolite, we propose a chemotactic mechanism of the fungal zoospores 
moving away from the metabolite gradients. There is strong experimental evidence 
of zoospore chemotaxis (see, e.g., Islam and Tahara [22]) and of negative zoospore 
chemotaxis for certain species (see [3] [7]). Moreover, it is well-documented that amino 
acids may act as chemo-repellent molecules for certain zoospores (see Nelson |34j). 
Under such considerations (even though it is not clear whether this biocontrol mech- 
anism is of chemotactic nature) we find it reasonable to assume that the inhibition 
process has a chemotactic component. In addition, up to the authors' knowledge, 
there are no reliable measurements of diffusion or chemotactic coefficients for the F. 
oxysporum in the literature. Thus, the purpose of the model is just to provide a sim- 
ple mathematical description of the underlying dynamics, which help us to explain 
qualitatively the emergence of stable invasive fronts for one species and to describe 
their effects on the biocontrol of the second species through very basic mechanisms. 

As a result, the model studied here consists of a three-component reaction-diffusion 
system of partial differential equations for the concentrations of the bacteria, the fun- 
gus, and the metabolite agent. It is assumed that the latter acts as a chemo-repellent. 
This fact is expressed through a chemotactic term of Keller-Segel type in the equa- 
tion for the fungus with the appropriate sign. In addition to diffusion and chemo- 
taxis, we incorporate logistic growth terms for both the bacteria and the fungus, as 
well as simple linear degradation and production terms for the chemical. Following 
well-established chemotactic cell-kinetic models (see p~9j[40]), we explore the small 
cell/zoospore diffusivity regime. This means that the diffusion coefficients of both the 
bacteria and the fungus are small relative to the diffusivity of the chemical. In par- 
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ticular, we found that the solutions reach stationary states in the absence of diffusion 
for the bacteria. These steady solutions approximate well the metastable patterns ob- 
served numerically when the diffusivity coefficient of the bacteria is taken to be very 
small. Although the metastable states eventually disappear and reach equilibrium 
solutions, they persist for very long times. (For instance, we conjecture that transient 
states resembling metastable patterns are what Swain and Ray [37] observed in their 
in vitro experiments.) It also to be noted that in both the zero and small bacterial 
diffusion regimes the chemotactic term prevents the invasion of moving fronts for the 
fungus concentration, simulating the biocontrol properties of the bacteria. Such fronts 
get stabilized by the interplay of diffusion, mean curvature and the negative chemo- 
taxis. We show asymptotically and numerically that, for a specific set of parameters, 
stationary or transient fronts emerge and reach equilibrium configurations. Thus, the 
simple mathematical model studied here captures the basic dynamical features of the 
antagonistic interaction of the two species observed in experiments. 

It is important to mention that, although the model was designed according to the 
experiment of Swain and Ray, it is applicable to any system of two antagonistic cell 
colonies, with one of them producing a chemo-repellent for the other. This paper joins 
the growing number of works studying multi-species chemotaxis models. Systems of 
Keller-Segel type with repulsion and attraction of multiple species have been mathe- 
matically analyzed by Horstmann in [21] ; the work by Conca et al. [14 describes the 
conditions for global existence and blow-up in the case of two-species and one chemo- 
attractant (see also the study of a one-dimensional attraction-repulsion system by 
Perthame et al. [36 ). One-species Keller-Segel systems with logistic growth terms 
have been studied by Mimura and Tsujikawa [30] and by Funaki et al. [18] . Winkler 
[4Tj 142] and Tello and Winkler [38] analyze the circumstances under which logistic 
growth can limit aggregation effects in one-species models with positive (attractive) 
chemotaxis. A two-species model (also of Keller-Segel type) with logistic growth has 
been recently analyzed by Kawaguchi [23], who studies a three-component system 
with a chemo-attractant for one of the species. It is to be observed, though, that 
Kawaguchi reports unstable two-dimensional front-, spot- and wave-like solutions in 
a chemo- attractive setting, producing branching instabilities. In contrast, this paper 
provides asymptotic calculations and numerical evidence of the stabilizing effect of 
chemo-repellent interactions, producing stable/metastable patterns which resemble 
the ones observed in certain experiments. The system studied here is, we believe, 
the simplest reasonable chemotactic model capturing the stabilizing effect on invasive 
fronts of the chemo-repellent interaction of two antagonistic species. 

The structure of the paper is as follows. In section [2] the model equations are 
presented. The derivation of the asymptotic equation of motion for an invasive front 
is given in section [3] The computation of stationary solutions in the zero-diffusion 
limit for the bacteria is the content of section[4] In section[5]we discuss the asymptotic 
conditions under which these stationary solutions are linearly stable. In section [6] we 
present the results of our numerical simulations and compare them to the asymptotics 
of the previous sections. Finally, section [7] contains a discussion of our results. 
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2 Modeling 



We shall denote by v the concentration of the bacteria which produces the chemo- 
repellent (metabolite) agent. The concentration of the latter will be denoted by c. 
The variable u denotes the concentration of the pathogenic colony (e.g. fungus). In 
order to describe the antagonistic activities (biocontrol) of these microflora colonies, 
we propose the following parabolic reaction-diffusion system of equations, 

u t = D u Au + Xu(uq — u){u — u*) — V • J c , 

v t = D v Av + /3v(v - v)(v - v*), (1) 

c t = D c Ac + Sv — ac, 

with (x,y,i) G Q x [0, +00), and where Q C M? is a bounded domain in the plane. 
According to custom A = d 2 + d 2 denotes the Laplace operator. The domain ft will 
be taken as a square or a circle with same area, that is 

n = [0, L) x [0, L) or n = {x 2 + y 2 < L 2 /tt}, 

being L > the characteristic length of the Petri dish. 

The concentration v has diffusivity equal to D v > and the value v = vq denotes 
the maximum sustainable population concentration of the bacteria. The produc- 
tion term is of bi-stable (or Nagumo) type [33, 29 , modeling logistic growth with a 
threshold. This enables the states v = and v = vo > to be stable equilibria in the 
absence of diffusion. We also assume that the unstable equilibrium point v* satisfies 
< v* < vo. Observe that the equation for v is decoupled from the other two. The 
chemical c is produced by the bacteria at the rate 6 > and it degrades at the rate 
a > 0. Its diffusivity coefficient is denoted by D c > 0. Finally, the equation for u 
is coupled to v and c via a chemotactic term J c . The concentration u diffuses with 
D u > and the reaction term is also of bi-stable type, with the roots slighted shifted 
from u = uq > and u = u* due to the presence of the chemotactic term, as we shall 
see later on. We assume, of course, that < u* < uq. 

The system of equations ([T]) is further endowed with no-flux boundary conditions 
of the form 

Vu • n = 0, \7v • n = 0, Vc • n = 0, at dtt, (2) 

with n being the unit normal to the boundary 90, reproducing the physical conditions 
of the experiments in vitro. 

The general form of the chemotactic flux J c is given by J c = — ux(c)\7c, where x 
is known as the chemotactic sensitivity function [31]. In this work we take 

X(c) = Xo, (3) 

with xo > constant, which is a customary choice when it is assumed that the 
species u always responds to a chemosensory stimulus in an uniform manner [3TJ [25] . 
Moreover, the parameter xo measures the intensity of the chemotaxis. The negative 
sign in J c indicates that the metabolite c acts, not as the customary chemo- at tract ant, 
but as a chemo-repellent agent blocking the spreading of the species u. This type of 
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movement toward lesser concentrations of a chemical is called negative chemotaxis. 
In the case of zoospores, such phenomenon is well-documented [3|, although it has 
not been as thoroughly studied as its bacterial counterpart. Therefore, we choose the 
simplest uniform chemotactic response function (|3|, which is also the first term in 
the expansion for small concentrations of the more general Lapidus- Schiller kinetic 
receptor law [26j[l7], commonly used for bacteria. Notice that, assuming ([3]), the 
chemotactic term in the equation for u in ([!]) has an advection component of the form 
XoV^x- Vc and an apparently less dominant kinetic-type term of the form xo^Ac. The 
effect of the latter is a slight shift on the equilibria when c reaches a stationary state. 
We observe that, although the system Q models the evolution of three components, it 
is essentially an scalar model for the fungus concentration u, because the equations for 
c and for v are independent of u. In this fashion we capture the simplest mechanism 
of the inhibition of the fungus due to the presence of the colony v producing a chemo- 
repellent c. 

Finally, we build the non-dimensional version of system ([!]) by making the following 
substitutions 
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D u D v I a . , 

D u — >• — — , D v — >• — — , t — » at, x \ — —x, — —y, 

CUV U* V* 

c— >> — , u — » — , v^ — , u* — ^ — , v* — » — , 
c u v u v 

a a coa ZJJ C 

The resulting non-dimensional system reads 

u t = D u Au + Xu(l - u)(u — u*) + xoV • (wVc), 

v t = D v Av + f3v(l - v)(v - v*), (4) 

c t = \ Ac + Sv — c, 

for (x,y,t) G x [0, +00), and where the unstable equilibria satisfy < u*,v* < 1. 
Without loss of generality we assume that L^/2D c /a = 0(1), so that the domain will 
be taken as 

n = [0,1] x [0,1], or n = {x 2 + y 2 < 1/tt}, 

namely, the unit square or the circle with area equal to one. 

In the sequel we analyze solutions to system Q. We are particularly interested 
in the dynamics of a moving front for the u variable. The asymptotics and stability 
of such a front, as well as its numerical computation, is the content of the following 
sections. 



3 Asymptotic equation of motion for the front 

When diffusion is small (D u « 0), localized patterns or layers in the u variable 
can be well approximated by interfaces. Here we derive the evolution equation of 
an interfacial curve when the concentration u is near the equilibrium configuration. 
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Such equation will be controlled by the concentrations u and c, and the location of 
the interface is coupled with v via c. We assume that the width of the interface is 
small. Let us denote 

E(£) = {(x,y) e ft : u(x,y,t) = u 2 } 

as the moving front, so that the outer and inner regions are defined by 

ftin = {(x,y) e ft : u(x,y,i) < u 2 }, 
ftout = {(x,y) e ft : u(x,y,t) > u 2 }. 

The equilibrium point u 2 will be determined later. We describe the motion of the 
interface in local curvilinear coordinates with components ((x,y,t) and r(x,y,i), 
normal and tangent to the interface, respectively, and normalized such that |V£| = 
|Vt| = 1. When diffusion is small, the dependence of u with respect to the tangent 
component is negligible and we approximate 

u(x,y,t) « u(£(x,y,t)). 

Under the assumption that the concentrations v and c have reached stationary (or 
quasi-stationary) states, we denote their values at the interface as and cs, re- 
spectively. Upon substitution into the equation for u in Q we obtain an ordinary 
differential equation for the solution u, given by 

(-s + D u k - XoVC • Vcz)v! = D u u" + \u(u - - u) + xoAc s ^, (5) 

where s = —d t ( is the normal velocity, and k = — A£ is the local curvature (see, e.g. 
[35]). The dynamics of the front is governed by the time-independent values of v and 
c at E. 



3.1 Interface equation of motion 

The total velocity along the normal direction has two main contributions: the velocity 
of the front in one dimension, which is determined by the nonlinear Nagumo term, 
and the chemotactic velocity. In order to obtain the former, note that if s and k 
were independent of ( (constant in time) then equation ([5| would resemble a one- 
dimensional front equation of the form 

— s\v! = D u u" + Xu(l — u){u — u*) + xoAcs^. (6) 

Since the value of Ac^ does not depend on time (it acts as a coefficient), the reaction 
term of last equation is also of Nagumo type, resulting in a shift of the roots for the 
stable/unstable equilibria. Rewriting the right hand side of (J6|, we obtain 

— s\u = D u u" + \u{u\ — u)(u — u 2 ), (7) 

where the equilibrium points U\ and u 2 now depend on the value of Ac^ as follows, 

U! = \{1 + u*) + ^V(l-^) 2 +4 Xo Ac E /A^ 1 + W1 X ° . Acs, 

/ Xo { ' 

u 2 = ^(1 + u*) - ^ V (1 - ii*) 2 + 4xoAc s /A « ^* - — rAc s . 

A^l — iz* J 
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Although the velocity of a one-dimensional front for a general reaction function 
can be obtained explicitly via a variational characterization [5], the Nagumo speed 
for this modified reaction term is directly computable. Indeed, equation ^ is an 
ordinary differential equation for a monotonic profile u(-) connecting the two stable 
equilibria u = and u = u\ . Suppose that u connects these equilibrium points on the 
left, that is u(— oo) = u\ and iZ(+oo) = 0. Therefore the profile is monotonic with 
v! < 0. In the phase space the solution to ^ is given by u' = (j){u), which leads to 
an equation for 0, namely 

—si(j> = Ax^-rr + Xu(ui - u)(u - u 2 ), 
au 

with boundary conditions (f)(0) = <j)(u\) = 0, and subject to the constraint <p < 0. 
The solution is of the form <j)(u) = —jiu(ui — u) with \i > 0. Since the velocity is 
independent of </>, upon substitution we obtain \i — y/X/2D u , yielding in turn 

si = y/2\D u (± Ul - u 2 ). (9) 

The sign of s± is that of u\ — 2u 2 ~ 1 — 2u* + 3xoAc$]/A(l — u*) in view of ([8|. 

Continued inspection of equation ^ shows that there is a contribution to the 
speed of the front due to the chemotaxis. Such term is known as the chemotactic 
velocity [TTJ [24] . In our two dimensional setting it has the form 

_ . _ dc 
d( 

In this fashion, we obtain the interface equation of motion (similar to that of [35] , 
but with a chemotactic contribution): 

dc 

s = si ~Xo-^ +D u k. (10) 

Remark. We observe that in the absence of chemotaxis the front will stop whenever 
Si = 0. The one dimensional case for a single Nagumo equation was studied by Chen 
et al. [13] when s\ was taken to be a function of the space variable. It was shown 
that when s[ < at the steady state x for which s(x) =0 the front is stable. In the 
present case we do not consider the two dimensional analogue of this situation which 
arises when the c variables acts as a death factor for the u equation, increasing the 
death rate of the pathogenic fungus. Instead, we concentrate on the new role of c as 
a chemorepellent. 

3.2 The case of a circular front 

Let us examine the case of a circular front, for which the interfacial curve can be 
written as 



£(*) = {(£, y) e Q : u(x, y, t) = u 2 } = {{x, y) e : V% 2 + V 2 = R(t)}, 
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where the function R : [0, +00) — >• [0, +00) indicates the localization and time- 
evolution of the front. In polar coordinates with r — y / ' x 2 + y 2 , the normal com- 
ponent is £ = R(t) — r , so that the curvature is n = — A£ = 1/r and the normal 
velocity is s = —d t ( = —R(t). Substituting into (10), the dynamics of a circular front 
is governed by the equation 

R(t) = -y/2XD:^u 1 - U2 )-(xo^ + — ) , (11) 

yz J \ dr r J \ r =R{t) 

which conveys the combined effects due to curvature, chemotaxis and the Nagumo 
speed. Here the value of s\ is the one computed in (|9|, because we are assuming that 
u connects u = u\ at r = +00, with u = at r = 0, in as much as ( points out in the 
direction of the origin and, consequently, dc/d( — —dc/dr. 



4 Stationary solutions 

In this section we assume that the concentrations v and c have reached stationary 
(in the case of zero diffusion D v — 0) or metastable (in the case of small diffusion 
< D v <C 1) states, which, in addition, have radial symmetry. For simplicity, we 
suppose that the spatial domain is the circle £1 = {x 2 + y 2 < 1/tt}. Since the equation 
for v is independent from the other two, let us denote by Voo(r) the radially-symmetric 
stationary solution to the equation for v in Q. Therefore, looking for a stationary 
solution for the chemical concentration c reduces to solving 

|Ac + ^oo-c = 0, in ft, 

(12) 

Vc • h = 0, at dfl. 
We first examine the case of zero diffusion for v. 

4.1 The zero diffusion limit 

In the limit of zero diffusion for the bacteria (D v = 0), the equation for v in Q is an 
ordinary differential equation whose solution tends to the stable equilibrium points 
v = or v = 1, depending on the initial spatial distribution. Let us suppose, for 
example, that the initial condition is a Gaussian of the form 

v(x,y,0)=Ae-" r \ (13) 



where r = ^ x 2 + y 2 and with A > 0, uj > 0. The mass of the initial condition shall 
be bounded below and above in order to allow the formation of a plateau inside the 
domain. Hence, we shall further assume that 

1< — < e w/ \ (14) 

Therefore, the stationary solution for v is determined by the plateau 

"-<•-> = (o bT.V < 15 » 

0, Ri < r < 1/V7T, 
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where the radius R\ satisifes < R\ < 1/ and is given by 



Ri 



]n(A/v*). 



(16) 



Then, the solution to (12) has the general form 
c(r) = 



C\I (V2r) + C 2 K (V2r), R 1 < r < l/^i, 
C 3 I (V2r) + 5, 0<r<R u 



(17) 



where 7 n , K n denote the modified Bessel functions for each n = 0, 1, . . . (see [HE]), 
and Ci, i = 1,2,3, are constants. Recall that Kq diverges at zero. The solution 
is subject to the Neumann boundary condition at r = and to ^-matching 

conditions at r = namely 



dr 



(dloiVtr) + C 2 K (V2r)) =0 



(CiJ (\/2r) + C 2 K (V2r) - C 3 I (V2r)^j 
^ (d I (V2r) + C 2 Jfo(>/2r) - C 3 / (v / 2- 



\r=Ri 



\r=Ri 



0. 



Using the known relations (see, e.g., [lj Iq(x) = h(x), Kq(x) = —Ki(x) and 
i^ n (x)/ n+ i(x) + K n +i(x)I n (x) — 1/x, we solve for d to obtain 



C 2 = y/25R 1 I 1 (y/2R 1 ), 
V25R! 



(18a) 
(18b) 



C 3 = r , * (gi(v^A)Ji(V2fli) - /i(v^)gi(v^fli)), (18c) 



completing the form of the solution (17). 



4.2 Small diffusion: metastable patterns 

When D v > 0, the solutions to the equation for v in Q with generic initial data and 
Neumann boundary conditions converge to the stable equilibrium solutions v = 1 or 
v = as t — )• +oo (see [lOj HI]). When diffusion is small (0 < D v <C 1), however, 
this convergence is very slow and the solutions can exhibit dynamic metastability 
[H |9]. After a short transition time dominated by the Nagumo term, the solu- 
tion v forms a pattern of layers which is apparently stable. These slowly evolving 
metastable solutions are neither local minimizers of the energy nor necessarily close 
to an unstable equilibrium solution. Although these solutions eventually decay to 
the patternless equilibria when t —> +oo, the time scale for substantial motion of 
the layers to occur increases exponentially depending on the size of the domain, the 
diffusion coefficient, and the size of the Nagumo term. For example, in the case of 
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the bistable react ion- diffusion equation in one dimension, namely u t = Du xx + f{u), 
Carr and Pego [8] estimated this mean transient time as ra(T) ~ exp(Cl/^D), where 
C = mm{f"(0), /"(l)}, and / > is of order of the minimum distance between layers 
and the distance between the layers and the boundary. In our setting, C = 0(f3) 
and / = 0(1), meaning m(T) ~ 0(exp(/3/^/D)). Even in one dimension, numerical 
estimations of the transient time seem to depend on D, the size of the domain, and 
the discretization mesh [T5] . 

Up to the authors' knowledge, there is no analytical estimation available for the 
mean duration of transient patterns for a bistable reaction-diffusion equation in a 
two-dimensional spatial domain^] Although we cannot provide an estimate for the 
life span of metastable solutions, our numerical simulations show the existence of a 
transient state in the v variable that induces, in turn, metastable patterns for the 
solutions of c and u. The stationary solutions for v and c of the previous section 
in the zero diffusion limit approximate well these metastable patterns (we refer the 
reader to section [6] for details). 



4.3 Approximate equilibrium solutions 

We discuss the conditions for the emergence of an equilibrium circular front in the 
variable u and provide an asymptotic formula to track its location. In this section we 
suppose that the diffusion coefficient for the species u is small, D u <C 1, so that the 



interface equation of motion (11) is valid in the geometric front propagation limit. 
As before, the diffussivity coefficient D v is either zero or very small, D v <C 1 (no 
relation between D v and D u is assumed at this point). Let us assume that r = Ro is 
the equilibrium interface position, with < Ro < 1/tt. Here we suppose that either 
R(t) — >• Ro when t — >• +00 (stationary case in the zero-diffusion limit for the bacteria, 
D v = 0), or that there exist To = 0(1/ X) > 0, T\ ^> To, and an uniform e > 0, such 
that 

\R(t)-R \ <e, for te [T ,Ti], 
in the metastable case when diffusion is small < D v <C 1. In both situations, 



the interface equation for motion (11) implies the following necessary condition for 
equilibrium, 

R(t)\r= Ro = -^2\D u (\ Ul -u 2 )-^- xoc'(Ro) = 0. (19) 

Ho 

In view that for small times metastable patterns resemble the stationary states in 
the zero-diffusion limit, we shall use the stationary solutions ( fl5| and ( [l7| ). Therefore, 
the stationary value of Ac^ at equilibrium is given by 2c(Rq) — Sv QO (Ro). After 



making uj large enough in the initial distribution (13) for v, we assume, without loss 
of generality, that the equilibrium radius is outside the plateau, i.e. that R\ < Ro. 
Thus, v OQ (Ro) = and we approximate 

Acx\ r =R « 2cOR ). (20) 



1 Numerical estimations for a square domain, however, have been provided by Horikawa |20| in 
the case D = 1, showing that log 10 m(T) ~ 0(g(l)), where g = g(l) is a linear function of I > 0, and 
/ is of order the size of the domain (see figure 6 in |20|V 
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Hence, we approximate the shifted equilibrium points u\ and ii 2 at the equilibrium 
position using (J8|, that yields 

Ul ~ 1 + V c (^o), ^2 ~ ^* - ttt: rc(i?o). 

A(l - ix*) A(l — u*) 

From ((TrJ) and ([l8|, we have that c(R ) = dIo(V2Ro) + ^^(V^^o), and 
c f (R ) = V2(Ci/i(v2Ro) - ^^(V^^o))- In view that # is bounded below, we 
use the approximations [4] Io(x) « 1 + x 2 /4, « x/2, if o 0*0 ~ ln(2/x) — e and 

i^i(x) « 1/x (where e « 0.5772 is the Euler constant), in order to estimate 

c(iJo) + C 2 (ln V"2 - e - In # ), 

c'(Rq) ~ Cii?o — C2/R0' 



Therefore, 

- (Ci (1 + iiZg) + C 2 (ln V2 - i - In R )) . 



1 1 3xo 

U 2 - TyUl « 1£* 



2 A(l-u*) 



Substituting into (19), and multiplying by i?o, we obtain 

p(Ro) = 0, (21) 

where the function p(-) is defined by 

p(x) := a3^ 3 + a 2 £ 2 + ai# + ao + foe In 

with coefficients 

a3 "V2A(1-^) >/2A(l-ix*)' 

«2 = X0C1, 

ax = ^2AA:(i - ^ + W1 3X ° r (d + C 2 (ln V^2 - 
V2 A(l — 

ao = A* - XoC 2 , 
= -y z\D u — = — -. 

A(l-U*) y/\(l-u*) 



Equation (21) provides an asymptotic formula for the equilibrium position. One 
may apply Newton-Raphson method to compute the zeroes of p in the interval 
[0, Once the approximated equilibrium position Rq has been computed nu- 
merically, one may use the asymptotic expansions I(J (x) ~ ^(l + 3x 2 /8) and Kq(x) « 
1/x 2 , to estimate c"(R ) « Ci(l + 3ii§/4) + C 2 /i?g. This yields, 

& - Xoc"(i?o) « 4 (A, - XoC 2 ) - X0 Ci(l + \Rl). (22) 
it /t 



The sign of the expression in (22) plays an important role in the stability of the front 
as we shall see in the next section. 
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We have thus constructed a steady state solution where the concentration for the 
variable u equals the equilibrium u = u\ « 1 outside a circular region of radius 
Ro determined by the solutions of the geometric front propagation equation. The 
chemical c is concentrated inside the plateau of radius R\ and its gradient balances 
the invading velocities. We will see how this basic state and its stability can explain 
the behavior of the solutions under biologically relevant initial conditions. 



5 Linearized stability 



To study the stability of the invasive front we consider again equation (19). Once 
again we suppose that both D u and D v are small, and D v possibly zero, with no 
relation between the two. The steady state solution outlined in the previous section 
can be perturbed in the form Rq + 7j(i), with rj(t) <C Ro- Linearization of equation 



([19]) gives 

(Du 



v(t)=( 1 ^-Xoc"(Ro))r ] - (23) 

This equation shows that the equilibrium point is linearly stable relative to radial 
perturbations of the front provided that 

X0 c"(R ) > ^. (24) 

If c"{Rq) now becomes a slowly evolving function of time, then it follows that, as 
c"(Ro) diminishes, the equilibrium point will lose stability and the front will start 
propagating. This will be the case for slowly diffusing chemo-repellent bacteria. As 
they diffuse with < D v <C 1 the concentration plateau will expand and the equation 
for c will have a flatter solution as the plateau moves. This will cause a smaller 
concentration gradient until the final steady state is homogeneous. This prediction 
is, in fact, verified by the numerical solution which will be discussed in the following 
section. 

The next question is the stability relative to azimuthally dependent perturbations. 
Again we consider the limit of a moving front and use the equation of motion driven 
by mean curvature for the interface. If the interface is parametrized by X(6,t) = 
(x(9,t),y(9,t)), where is the polar angle, then the equation of motion for the front 
takes the form [16] 

X t = (c- D u K(0,t) -Vc(X(0,t)) • n)n, (25) 
where n is the outward unit normal vector, and 

m x _ xeyee Ve^ee 

K{ , ) ~ (4 + yl) 3/2 

is the local curvature. The steady states are given by 

c - D u k(6) - Vc(X(<9)) • n = 0, (26) 
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which provides an ordinary differential equation for the curve X(0). In the present 
situation it is equation ( [19] ) . To consider the full linearized stability problem, we 
take X(0,i) = Xq(0) + £(0,£), where Xq{0) is the solution of the steady state. The 
linearized equation for the perturbation takes the form 

( 3Ko II " ' ^ ~ X ° ' ( ^)™ " Xoie > (^H), (27) 



where no is the outer normal to the curve Xq(0), the ' denotes derivative with respect 
to 0, and (Xq) 1 - and (£ / )~ L are orthogonal to Xq and respectively. The matrix D 2 c 



is the Hessian of the concentration at the equilibrium front. Equation (27) provides 
the parabolic system for the perturbation front £. These equations can be further 
simplified by choosing a representation for £ in the coordinate system formed by the 
unit tangent fo and the unit normal no to the unperturbed front Xq. We thus take 

Z = p(0,t)To(0) + q(0,t)no(0). 



It follows immediately from (27) that p t = and hence we may take p — 0, in as 
much as perturbations in the tangential direction amount to a change in parametriza- 
tion. We are thus left in the present case with an equation for q in the form 

qt = ^qee + (Jr ~ Xoc"(R ))q, (28) 



with periodic boundary conditions in 9. Equation (28) shows that azimuthally per- 
turbations decay if the linearized stability condition (24) also holds. It is to be noted 
that, for an arbitrary azimuthally dependent steady state, the stability equation will 
be a scalar equation with ^-dependent coefficients. The stability will be determined 
by the spectrum of the associated differential equation. In particular, ( [28] ) shows that 
fronts with sign changing chemotaxis may be stabilized by diffusion. On the other 



hand, equation (28) may be helpful to explain instability results (see [23 ) for a related 
system with the opposite chemotaxis sign. 

The issue of nonlinear stability will be addressed numerically in the following 
section. 



6 Numerical results and nonlinear evolution 

To study the nonlinear evolution of solutions we solve system Q numerically. For that 
purpose, we implemented an explicit finite-difference Euler scheme in an appropriately 
scaled square domain 0<#<1,0<2/<1, using a grid spacing with Ax = Ay = 
« 3.9 x 10" 3 and N = 256. The time step was set as At = 10" 3 x Ax w 
3.9 x 10 -6 . Hence, the Courant number /i = At/ (Ax) 2 is approximately fJL ~ \, 
assuring numerical stability. In addition, we imposed zero-flux boundary conditions 
for all species. 

Although explicit Euler schemes are not commonly used to solve stiff equations 
(because they require very small time steps in order to avoid instabilities), they are 
particularly easily implemented to run in parallel on a Graphics Processing Unit 
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Table 1: Non-dimensional parameter values used in the computations. 



Description 


Symbol 


Value 


Diffusion coefficient of u 


D u 


0.01 


Diffusion coefficient of v 


D v 


10~ 5 or 


Rate of production of c 


S 


10.0 


Reaction coefficient of v 


P 


8.0 


Reaction coefficient of u 


A 


60.0 


Chemotactic sensitivity 


Xo 


3.2 


Unstable equilibrium for u 


u* 


0.2 


Unstable equilibrium for v 




0.5 



(GPU). In this way, the drawback of a small time step is compensated by the high 
performance parallel computation capabilities of a GPU with hundreds of processors. 
Our computations were performed on a commodity-type NVIDIA® GeForce GTX 480 
graphics card with 480 stream processors. In this fashion, we were able to compute 
tenths of millions of time steps in a few hours. 

Finally, to perform the numerical calculations it is necessary to determine the 
non-dimensional parameter values involved in model Q . Table [I] shows the set of 
parameter values used during our simulations. 



6.1 Zero-diffusion 

We begin by studying system Q with D v = 0. The initial conditions for the bacteria 
(species v) is a colony concentrated at the origin of the domain, corresponding to a 



centered Gaussian of the form (13), with A = 3, oo = 1000 and v* = 0.5, that is 



v(x, y, 0) = 3e -iooo((,-o.5) 2 + ( y -o.5) 2 ) _ (29) 

It produces a chemical c which initiates with concentration c(x, y, 0) = 0. The initial 
value for the species u is 

u(x,y,0) = i (>°°((-°-2) 2 +(y-o.2) 2 ) +e -iooo(( a; -o.2)^ + ( !/ -o.2)^)j- 1 ) (30) 

which represents a localized colony at one of the corners of the domain. 

The evolution of the solution to system Q is displayed in figures [l] to [3] In figures 
l(a)||l(b)| l(c)| and |l(d)| we observe the formation of a sharp plateau for the bacteria 



concentration v. The cross section of the plateau has radius given analytically by 



(16) and, with the parameter values A = 3, uj = 1000 and v* = 0.5, then it has 



an approximated value of R\ = 0.0423. (We can also evaluate the constants in 



(18) to obtain a theoretical stationary solution c given in (|17|)). In these numerical 



simulations, the computed value for the radius of the plateau Ri was approximately 
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Ri ~ 0.0412 , which gives a relative error of 2.6% with respect to its theoretical value. 
The chemical concentration c reaches a steady state which is displayed in figure l(e)| 
at time t = 9.8039. 



Figures 2(a) and 2(b) display the formation of the initial invasive front in the vari- 
able u which generates from the initial condition. Note that it becomes perpendicular 
to the boundary due to the Neumann conditions. At the time steps of figures 2(c)| and 
|2(d)| the plateau of bacteria and the steady state for the chemical are already formed, 
causing the emergence of a repelling region in the center of the domain for the front 
in the variable u. In figures 2(e) and 2(f) it is shown how the front responds to the 
chemical gradient by bending according to equation ( [25] ) . Finally, the front encircles 
the repelling region forming the steady state whose depression is centered around the 
bacteria. The steady state for u is depicted in figure |2(h)| at time t = 9.8039. All 
surfaces in figures [2(a)] - 12 (h)| are colored according to the chemical concentration for 
an easier interpretation of the chemotactic effects. 

Finally, figure [3] shows a top view of the final steady state for the front. Recall 
that we may compute a theoretical radius Ro of equilibrium by solving equation (21). 
With the parameter values depicted in Table [I] and taking A = 3, u = 1000, and 
S = 10 in (13), there is only one zero of p(-) in [0, 1/tt], and hence, the predicted 
equilibrium radius is approximately 



Ro « 0.1315. 



(31) 



It is to be noticed that the computed numerical value for the equilibrium radius 
of the solution depicted in figure [3] is Rq « 0.1235. A comparison with the theoretical 
value in (31) gives a relative error of 6%. Moreover, in both cases the stability 



condition (24) holds. For example, for the theoretical value of Rq we obtain 



Xoc"(R ) « 3.4409, 



IK 



0.5783, 



(32) 



yielding the negative sign of (22) in order to fulfill the stability condition. It is 



important to remark that condition (24) is necessary for stability of the front under 



both radial and azimuthal perturbations. 



6.2 Small diffusion: < D v < 1 

In order to explore numerically the effects of D v ^ 0, we took D v = 10 -5 as the 



diffusion coefficient for the bacteria, together with the same initial conditions (29) 



and (30) for the bacteria and the fungus, respectively. The initial distribution for 
the chemical was taken, once again, as c(x,y,0) = 0. In such a case, the bacteria 
will themselves evolve according to a spreading front in the bistable reaction-diffusion 
equation. It is known [lOj [IT] that under Neumann boundary conditions and positive 
diffusion, the solutions will eventually diffuse and decay to one of the two stable 
constant equilibrium solutions v = or v = 1. Under the parameter values of our 
simulations and the initial condition considered, we observed the case in which the 
bacteria concentration eventually vanishes for very large times as it is shown in figure 

m 
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After a short time, however, a metastable solution is formed. This pattern is 
smooth due to diffusion but resembles the numerical plateau of figure [T] Ob serve 
that this metastable state is formed at a time before t = 0.980 (figure 4(a)| ), and 
persists for a very long time of order t = 100 (figure 4(c)). The concentration of 
bacteria eventually vanishes (figure [4(d) [ ) and, before a time of order t = 180, it has 
already reached the steady state v = 0. 

As this behavior occurs, the concentration c of the chemical also forms a slowly 
evolving metastable pattern. Figure |5(a)| depicts this solution, which persists for a 
very long time. Notice the resemblance of the steady concentration c when D v = 



of figure 1(e) with the concentration when D v = 10 5 for small times shown in figure 
5(a)| After a long time of ordet t — 150, the chemical c reaches the steady state 



c = as it is shown in figures |5(b)| |5(c) and |5(d) In this case the chemotactic 



gradient vanishes and, as expected from equation (23), it is shown in figures [6(a)] to 



6(f) that the front of the fungus eventually covers the domain since the steady state 



is unstable. As D v becomes smaller, the time of existence of the steady (metastable) 
state is longer. 

It is important to remark that the computed solutions in the D v = case are good 
approximations of the metastable states when we perturb the diffusion parameter D v . 
Figures [7| [8] and [9] show the corresponding cross sections for the final steady states 
in the zero diffusion limit and the metastable states with D v = 10 -5 for short times. 
For instance, figure [7] shows cross sections of the concentration of the chemical c. 
The continuous (red) plot corresponds to the analytic steady solution (17) computed 
for a circular domain in the zero diffusion limit D v = 0. The dotted graph (blue) 
shows the numerically computed solution c on a square domain with D v = at 
time t = 9.8039, when the steady state has been already reached. Notice that it 
approximates well the analytic stationary solution. The error, which can be observed 
in the tails of both graphs, is due to the Neumann conditions computed on a circle 
(analytical solution) and on a square (numerical solution). The dashed (green) graph 
depicts the concentration c computed numerically with D v = 10 -5 at time t = 9.8039, 
when a metastable state has been reached, but in a short time step relatively to the 
relaxation time with small diffusion. It can be observed that the numerical solutions, 
both in the zero diffusion and small diffusion limits, approximate well the analytical 
solution (17). 

In the same fashion, figure [8] shows the cross sections for the concentration v of 
bacteria. The solid plot (red) shows the plateau ( 15 ) with an analytical radius equal to 
Ri = 0.0423. The dotted (blue) graph represents the numerically computed plateau 
with D v = at time t = 9.8039, with radius Ri = 0.0412. The dashed graph (green) 
shows the concentration for v at time t = 9.8039 when computed numerically with 
diffusion equal to D v = 10 -5 . It is smooth due to diffusion and it corresponds to a 
metastable state which remains for a long time, until it eventually reaches a constant 
state v = 0. 

Figure [9] shows the cross section for the numerically computed concentrations for 
u in the zero diffusion case D v = (continuous plot in blue), and in the case of small 
diffusion D v = 10 -5 (dashed plot in green). The former corresponds to a steady state, 
formed within a short time, and stable in the absence of diffusion for the bacteria. 
The latter corresponds to a metastable state for small diffusion. 



16 



6.3 The experiment of Swain and Ray 



Finally, we reproduce numerically the experimental observations of Swain and Ray 
[37] . To this end we take, as in the experiments, the initial concentration of the fungus 
to be uniform and constant, with value slightly above the threshold concentration 
u* = 0.2. Hence, our initial condition was taken as u(x,y, 0) = 0.21. Now we plant 
two localized colonies of bacteria in the form 



V(x,y,0) = 3 ( e -1000((x-0.2) 2 + (!y -0.5) 2 ) + e _1000((x-0.8) 2 + ( y -0.5) 2 )^ 



(33) 



simulating the administration of the bacteria B. subtilis into the control of the F. 
oxysporum. The initial concentration of the chemical was taken as c(x, 0) = 0, as 
before. In this numerical simulation we set the diffusion coefficient for the bacteria 
as D v = 0. 



Figures 10 and 11 show the dynamics of the solutions under these initial conditions. 
As expected, the bacteria and the concentration of the chemical evolve practically 
independently of the two colonies, as shown in figures [10(a)) [To(b)| [10(g) [To(dJ 



and 10(e) Notice the formation of two localized steady plateaus for the bacteria 
concentration, as shown in figure |10(d)| It induces a steady chemical concentration 
shown in figure [10(e) As in the previous example, the invasive front for the variable 
u encounters the chemical gradient, localized around two points of the domain. The 
repulsion causes again the front to go around the higher chemical concentrations 
depicted in figure 10(e)| settling into the superposition (to leading order) of two 
circular fronts in equilibrium, resembling the one studied in the first example. This 
evolution is shown in figures pT(a)] through [IT(h)| 

Finally, figure [12] shows the top view of the stationary state for the fungus. Notice 
the remarkable resemblance with the experimental pattern of Swain and Ray (figure 
2 in [37 ). The steady state computed here is a good approximation of the metastable 
pattern which occurs when diffusion of the bacteria is small. We conjecture that the 
latter is a suitable model for the experimental quasi- stationary state observed in vitro. 



7 Conclusions 

In this study we have shown how the process of inhibition of an invading front of 
one species, triggered by the chemical produced by another species, can be accurately 
described by two reaction-diffusion equations of Nagumo-type which are chemotacti- 
cally coupled to a third equation for the chemical. We have exhibited the existence of 
a radially symmetric steady state (for a circular domain), that was shown to be stable 
in the geometric front propagation limit. It was shown numerically that this steady 
state is an excellent approximation to the steady state obtained in a square domain 
as the result of the repulsion by the chemical of an invading front. The asymptotic 
solutions for the stability of the front suggested instability of the circular steady state 
as the chemical gradients decrease. This prediction is tested numerically with the 
expected results. 

In addition, we showed numerically that different colonies of bacteria act inde- 
pendently, producing a final fungus pattern which is well approximated by the su- 
perposition of the basic state studied here. This result explains qualitatively the 
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mechanism observed in the in vitro experiments performed by Swain and Ray [37] . 
We have thus shown how the solution of a system of equations can be understood 
in simple terms using the ideas of geometric front propagation, which were able to 
predict both quantitative and qualitative behaviors. Finally, we anticipate that a dis- 
tribution of attracting and repelling chemotaxis can produce dendrite-like patterns 
which can reach a steady state. 
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(a) t = 0.1961. 



(b) t = 0.9804. 




(e) Chemical concentration c at time t = 
9.8039. 



Figure 1: Figures (a) - (d) show the concentration v of the bacteria colony as time 
evolves. The bacteria v reaches a stationary state given by the plateau in figure (d). 
In this figures, the concentration for v is colored according to the concentration of 
the chemical it produces (see color chart on the right). The steady state for the 
concentration of the chemical c at time t = 9.8039 is shown in figure (e). 
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(a) t = 0.1961. 



(b) t = 0.7843. 




(g) t = 3.9216. (h) t = 9.8039. 



Figure 2: Concentration u of the harmful colony (fungus) for different times in the 
zero-diffusion limit for v (i.e. with D v = 0). Figure (a) shows the concentration of 
u at time t = 0.1961, near the initial condition p0| ). The latter was localized near 
one of the corners of the domain. Observe that u diffuses, senses the chemo-repellent 
localized in the center of the domain, and reaches a stationary state depicted in figure 

(h). 
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Figure 3: Stationary end state for the concentration u of the harmful colony (fungus), 
at time t = 9.8039. It is coloured based on the concentration c. 
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(a) t = 0.980. 



(b) t = 9.804. 




(c) t = 99.020. (d) t = 166.667. 



Figure 4: Figures (a) - (d) show the concentration v of the bacteria colony as time 
evolves when D v = 10 -5 . The concentration of v is now a smooth function (due 
to diffusion) that resembles the plateau of figure [T(d)| when D v = for short times. 
Figure (a) depicts this metastable state. The transient, however, tends to the equi- 
librium point v — when t — >• +oo. In this figures, the concentration for v is colored 
according to the concentration of the chemical it produces (see color chart on the 
right). 
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(a) t = 9.8039. 



(b) t = 49.020. 




(c) t = 88.235. (d) t = 147.159. 



Figure 5: Figures (a) - (d) show the concentration c of the chemical for different 
times, when D v = 10 -5 . Although eventually c reaches the constant steady state 
c = 0, figures (a) and (b) show the metastable state for short times. Note that figure 
(a) resembles the stationary state computed when D v = in figure 1(e) 
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(a) t = 0. 



(b) t = 1.961. 





(c) t = 4.902. 



(d) t = 5.882. 





(e) t = 50. 



(f) t = 99.020. 





147.059. 



(h) t = 166.667. 



Figure 6: Fungus concentration u for different times when D v = 10 -5 . Observe the 
emergence of a metastable pattern which resembles the steady state computed with 
D v = (figure [2(h) >. This metastable solution changes very little between times 
t = 5.882 (figure (d)) and t = 50 (figure (e)). Due to bacterial diffusion, the chemical 
concentration decays allowing the front to invade all the domain as shown in figures 
(f), (g) and (h). 
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Figure 7: The continuous (red) plot shows the cross-section for c computed analyt- 
ically in (17) in the zero diffusion limit D v = 0. The dotted plot (blue) shows the 
numerical solution for c when D v = at time t — 9.8039, when it has reached the sta- 
tionary state. The dashed line (green) shows the numerically computed concentration 
c when D v = 10 -5 , also at time t = 9.8039, depicting the metastable state. 



Figure 8: The continuous (red) plot shows the cross-section of the stationary state 
(plateau) v computed analytically in ( fl5| ) in the zero diffusion limit D v = 0. The 
dotted plot (blue) shows the numerically computed plateau for v when D v = at 
time t = 9.8039. The dashed line (green) represents the numerically computed con- 
centration v when D v = 10 -5 , also at time t = 9.8039 (metastable state). 
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Figure 9: The continuous (blue) plot shows the cross-section of the numerically com- 
puted stationary profile for the concentration u in the zero diffusion limit D v = 
at time t = 9.8039. The dotted plot (green) represents the numerically computed 
concentration u in a metastable state at time t = 9.8039 when the diffusion coeffi- 
cient < D v = 10 -5 is small. Observe that the stationary value computed in the 
zero diffusion limit is a good approximation of the metastable concentration for short 
times. 
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(e) Chemical concentration c at time t = 
9.8039. 



Figure 10: Concentration v of the bacteria colony for different times. It reaches a 
stationary state given by the two plateaus in figure (d), and induces, in turn, the 
steady state for the chemical concentration c shown in figure (e). 
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(a) t = 0. 



(b) t = 0.0392. 




(c) t = 0.1176. (d) t = 0.1961. 




(e) t = 0.2549. (f) t = 0.3137. 




(g) t = 0.4902. (h) t = 9.8039. 



Figure 11: Concentration u for different times when D v = 0. The initial condition 
shown in figure (a) is the uniform concentration u = 0.21. The fungus diffuses and 
eventually reaches a stationary state shown in figure (h). Notice that as it evolves, 
it avoids the repelling region induced by the chemical gradient produced by v, whose 
concentration is shown in figure |10(e)| The result is the formation of two equilibrium 
fronts. 
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Figure 12: Top view of the stationary end state for the concentration u of the harmful 



colony, at time t = 9.8039 (see figure [11(h) ). It is coloured based on the concentration 



c. Notice the remarkable resemblance with the experimental result of Swain and Ray 
(figure 2 in 
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